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ICASE Fluid Mechanics 


Due to increasing research being conducted at ICASE in the field of fluid mechanics, 
future ICASE reports in this area of research will be printed with a green cover. Applied 
and numerical mathematics reports will have the familiar blue cover, while computer science 
reports will have yellow covers. In all other aspects the reports will remain the same; in 
particular, they will continue to be submitted to the appropriate journals or conferences for 
formal publication. 
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Nonlinear stability of oscillatory core- annular flow: A generalized 
Kuramoto-Sivashinsky equation with time periodic coefficients. 

Adrian V. Coward * Demetrios T. Papageorgiou * Yiorgos S. Smyrlis * 


Abstract 

In this paper the nonlinear stability of two-phase core-annular flow in a pipe is examined 
when the acting pressure gradient is modulated by time harmonic oscillations and viscosity 
stratification and interfacial tension is present. An exact solution of the Navier-Stokes equa- 
tions is used as the background state to develop an asymptotic theory valid for thin annular 
layers, which leads to a novel nonlinear evolution describing the spatio-temporal evolution of 
the interface. The evolution equation is an extension of the equation found for constant pressure 
gradients and generalizes the Kuramoto-Sivashinsky equation with dispersive effects found by 
Papageorgiou, Maldarelli & Rumschitzki, Phys. Fluids A 2(3), 1990, pp. 340-352, to a similar 
system with time periodic coefficients. The distinct regimes of slow and moderate flow are con- 
sidered and the corresponding evolution is derived. Certain solutions are described analytically 
in the neighborhood of the first bifurcation point by use of multiple scales asymptotics. Exten- 
sive numerical experiments, using dynamical systems ideas, are carried out in order to evaluate 
the effect of the oscillatory pressure gradient on the solutions in the presence of a constant 
pressure gradient. 
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1 Introduction 


Core annular flows (CAFs) are immiscible two-fluid flows in a cylindrical tube in which one fluid 
moves through the tube core and the other liquid occupies the annular region surrounding the core. 
If a constant pressure gradient causes the flow, an exact solution of the Navier-Stokes equations 
is available with the interface between the two fluids being perfectly cylindrical. A considerable 
amount of research has been applied in an effort to characterize the instability of such perfect CAF 
arrangements. There are many technological applications where the CAF serves as a useful model 
for the dynamics, such as, lubricated pipelining, concurrent flows in packed beds, coating processes, 
and liquid-liquid displacements in the presence of a wall-wetting layer in porous media, to mention a 
few. The important theoretical question lies in the prediction of whether the core-annular geometry 
will be realized or whether the flow will break up into a slug or emulsion flow. 

Experimental and theoretical studies indicate that the main physical parameters that affect 
the stability are capillary forces, viscosity differences and density differences between the phases. 
The effect of viscosity differences was first studied by Yih (1967) for the stability of two-phase 
plane Couette-Poiseuille flow; he derived analytical expressions for the growth-rates of linear per- 
turbations in the long-wave limit (i.e. disturbances with wavelengths which are much larger than 
the plate separation distance) from which the following conclusions are inferred: If the two fluids 
have different viscosities an instability is possible at any value of the Reynolds number, however 
small (the instability is not associated with high Reynolds numbers alone). The flow is stable (for 
long waves at least) if the less viscous fluid occupies the thinner of the two layers and unstable 
for converse arrangements. Hooper (1985) and Hooper k Boyd (1983) generalized Yih’s results for 
semi-infinite and infinite two-phase Couette flows which are useful models in studying short-wave 
instabilities. Renardy (1985) addressed the problem numerically and shows that in planar Couette 
flow short and intermediate waves can become unstable in parameter regimes where the long wave 
analysis predicts stability. 

In cylindrical geometries the situation is different in that surface tension acts to destabilize 
interfacial waves longer than the core diameter even in the absence of flow (Tomotika (1935), 
Goren (1962)). For flow in a core-annular arrangement, Hickox (1971) analyzed long wavelength 
perturbations and found instability whenever the annular fluid is more viscous than the core fluid. 
More complete linear stability analyses have been carried out in a series of articles by Joseph, 
Renardy k Renardy (1984), Preziosi, Chen k Joseph (1989), Hu k Joseph (1989), Hu, Lundgren 
k Joseph (1990), Chen, Bai k Joseph (1990), Bai, Chen k Joseph (1992), Chen k Joseph (1991) k 
Chen (1992). These authors use numerical and long-wave techniques to produce a detailed picture 
of the linear stability to mainly axially symmetric perturbations. The authors find that there are 
essentially two competing mechanisms at play: (i) Capillary forces which destabilize long waves 
and are, to leading order in the interfacial deflection, independent of the basic flow, (ii) Viscosity 
differences which provide a jump in the leading order axial velocity perturbation which destabilizes 
short waves and can also destabilize long waves when the annular fluid is the more viscous (long 
waves are stabilized by these forces if the core fluid is more viscous). The studies detailed above 
show that the interaction of these mechanisms can produce a window of linear stability which can 
be explained as follows: in situations when the annular fluid thickness is much smaller than the 
core radius and its viscosity smaller than that of the core fluid, the viscosity difference mechanism 
provided by the axial velocity jump (see above) causes a stabilization of long waves (which would 
have been unstable otherwise due to the capillary instability) at sufficiently large Reynolds numbers. 
As the Reynolds number increases further the flow becomes unstable to the short wave mechanism. 
It is established, therefore, that the flow parameters can be chosen so that the flow is linearly 
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stable to perturbations of all wavelengths. For a complete and up to date account of theoretical 
and experimental results and their comparisons see the texts by Joseph & Renardy (1993). 

The nonlinear stability of CAFs when capillary forces are important, has been considered by 
Papageorgiou, Maldarelli & Rumschitzki (1990) (hereafter refered to as PMR). This was done in 
the limit of a thin annular layer by a systematic asymptotic expansion procedure on the Navier- 
Stokes equations. In this limit the film dynamics reduce to those given by lubrication theory and 
a matching can be achieved analytically between film and core dynamics to produce an evolution 
equation for the scaled interfacial amplitude 7]{t , z ) where t is time and z axial distance, which can 
be written as 


Vi + VVz + Vzz + VVzzzz + (1 - —)£(«) = 


( 1 ) 


where 0 < v < 1 is a parameter inversely proportional to the square of the wavelength and m is the 
viscosity ratio of core and film fluids. The last term arises due to viscosity differences and is a linear 
pseudo-differential operator with a known spectrum (see PMR and also below). The validity of the 
lubrication approximation is given by Georgiou, Maldarelli, Papageorgiou & Rumschitzki (1992) 
who analyze the linear stability dispersion relation in the limit of a thin annular fluid and identify 
the linear terms of (1) exactly to leading order. Further more, it is shown that the lubrication 
approximation agrees with the numerically computed eigenvalues for annulus to core radius ratios 
(denoted by e) as large as 0.2, a finding which implies that (1) may be a useful approximation for 
finite but small e (note that e can easily be re-introduced in (1)). 

When m = 1, (1) reduces to the Kuramoto-Sivashinsky equation which has been widely stud- 
ied both analytically (see for example Nicolaenko, Scheurer k Temam (1985), Constantin, Foias, 
Nicolaenko k Temam (1988)) and computationally (Sivashinsky k Michelson (1980), Hyman k 
Nicolaenko (1986), Hyman, Nicolaenko k Zaleski (1986), Kevrekidis, Nicolaenko k Scovel (1990), 
Papageorgiou k Smyrlis (1990) (referred to as PS), Smyrlis k Papageorgiou (1991) (referred to 
as SP) among others.) These studies show that the dynamics exhibits a low-modal behavior and 
complexity sets in as v decreases below unity. PS and SP showed numerically the existence of 
a period-doubling route to chaos according to the Feigenbaum scenario (see Feigenbaum (1980)) 
and have computed a Feigenbaum number from their numerical results with a three digit accuracy. 
It is noted that such results are reproducible with a fairly low- dimensional representation of the 
PDE; the number of modes required to capture the dynamics was found (empirically by extensive 
numerical experiments) to be a few more than , the number of linearly unstable 27T-periodic 
waves. 

An extension of the analysis of PMR to include asymmetric interfacial deflections and in par- 
ticular the effects of pipe rotation, has been carried out by Coward k Hall (1993). A spatially 
two-dimensional version of (1) is derived, analyzed and solved numerically. In the absence of ro- 
tation it is found that asymmetric initial conditions lead, after a long time, to axially symmetric 
motions. When rotation is present, however, the most unstable linear wave has azimuthal depen- 
dence and as a consequence the nonlinear solutions are two-dimensional. 

If the annular and core fluids have different viscosities (m ^ 1) the following behavior has 
been established in PMR for (1) and by Coward and Hall for the higher dimensional equation: 
Viscosity differences cause dispersive effects. If |1 — l/ra| exceeds a certain value (which depends 
on v), the long-time evolution produces nonlinear traveling waves whose amplitude grows with 
the amount of dispersion. In particular this is established for situations where in the absence of 
dispersive effects (m = 1) the KS equation produces spatio-temporal chaos. The conclusion is that 
dispersion (viscosity differences) can organize chaotic motions into well-defined stable travelling 
wave solutions. 
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The present work extends the research cited above to CAFs which have driving axial pressure 
gradients with time periodic modulations. Such flows are easily achieved in practice by using a 
pump to drive the flow. Time periodic two-phase flows have received very little attention. The 
stability of a single layer of fluid on a flat plate which is performing time periodic oscillations 
horizontally has been studied by Yih (1968) using a long wave analysis coupled with a Floquet 
theory. Von Kerczek (1987) studied the corresponding problem when the plate is vertical and 
performs oscillations along its plane, by solving the Floquet eigenvalue problem numerically for 
different wavenumbers. Both investigations indicate that the flow can be stabilized by the periodic 
oscillations. Coward &; Papageorgiou (1994) study the stability of two-phase Couette flow when 
the upper plate has a horizontal velocity which comprises of a constant component and a time 
periodic modulation. The problem is more involved due to the presence of two fluids, and analytical 
expressions have been found for the Floquet exponents (growth rates) in the limit of long waves. 
The main finding is that the oscillations can completely stabilize otherwise unstable flows, but at 
the same time can make unstable flows, in certain parameter regimes, even more unstable. 

This article considers the nonlinear stability of an oscillatory CAF. An exact laminar parallel 
flow solution of the Navier-Stokes equations can be found for an annulus of arbitrary thickness 
which bounds the core with a perfectly cylindrical interface. This time-periodic flow is used as 
the undisturbed state in the development of an asymptotic theory using e as the small parameter 
in order to develop a nonlinear evolution equation for the interfacial deflections which takes into 
account the effects of both capillary forces and the background harmonic fluctuations. The scaled 
evolution equation, using the notation of (1), which is valid for low frequency modulations, takes 
the form 

7} t + (1 + A COs(Ot))f??k + Vzz + V1)zzzz + (1 )£(«) = ( 2 ) 

Ttl 

where A is the nondimensional amplitude of the background pressure gradient modulations and 
0 their scaled frequency. This equation is studied analytically and computationally in order to 
characterize the effect of the background modulations on the dynamics. 

The paper is organized as follows. In Section 2 we formulate the problem, write down the 
unperturbed time periodic flow and give the exact nonlinear equations of motion along with the 
interfacial conditions. In Section 3 an asymptotic analysis of these equations is carried out in the 
limit e — * 0 and two canonical regimes are analyzed to produce equations of the form (2) with 
different pseudo-differential operators which correspond to slow flow with moderate surface tension 
and moderate flow with large surface tension, respectively. In Section 4 the equations are analyzed 
near the point v - 1 which corresponds to the minimum wavelength which first becomes unstable 
(the first bifurcation). A multiple scales analysis is employed to describe the solution and the 
analytical results are compared to numerical solutions with good agreement. Section 5 is devoted 
to extensive numerical experiments of (2) for values of v < 1 where more than just two or three 
modes are relevent. We address questions such as the effect of the background oscillations on (i) 
steady states of KS, (ii) time periodic states, (iii) chaotic states, (iv) travelling wave solutions, 
and (v) travelling nonlinear dispersive states. Section 6 contains the conclusions of this study and 
discusses future work. 

2 Formulation of the problem 

A fixed circular cylinder of constant radius r* - R 2 contains two incompressible, immiscible fluids 
which occupy concentric core and annular regions. These fluids have a common density /?, but 
different viscosities fi x and ^ 2 , where the subscript corresponds to regions 0 < r* < R\ and 
Ri < r* < R 2 respectively. 
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Polar coordinates are used so that r* = 0 is aligned parallel to the z axis and an unsteady axial 
pressure gradient 

VP* = - (F + A cos (u>t*)) z, 

is applied to both fluids. In the analysis that follows gravity is neglected if the pipe is horizontal 
and can be included quite easily in vertical arrangements (see [[ 8 ]], and [[ 1 ]]). 

The undisturbed state is a parallel time dependent flow of the form U 12 (r*,t*) = (o? 0? W 1)2 (^> r )) 
and is governed by following form of the Navier-Stokes equations in the core and annular layers 
respectively 


dW\ 

’ dt* 

dW* 2 

) 

dt* 


F + A cos (wP) + yui 
F + A cos ( u>t * ) + 


( d 2 W\ 1 dW[ 

f d 2 W 

l dr 


1 dW, 


2 , 1 2 

.2 ' 




The annular fluid satisfies a no-slip velocity condition at the cylinder wall and the core velocity 
must remain finite at the axis of the cylinder. The interface conditions dictate that velocity and 
tangential stress be continuous, so that 


Pi 


W m 1 {r m = R l ) = (r* = Pi) , 

dWl , * „ , dW* 2 , * _ , 

i(r* = B,) = Ii 2 ~zzr(r ‘ = Ri)- 


9r* v - "" “ ^ dr > 

The surface tension a between the fluids induces a jump in the normal stress at the interface, 

a 


p; = p; + 


Pi’ 


at r* = R\. 


These equations and boundary conditions admit the following solutions 

w, 

4^1 4 /j 

P j f . k 

[Bio {fhr‘) + CKo (Jhr") - l]exp (iwi*) f . 


* = F ... ’ ll + F + j— [AI 0 (ftr*) - 1] exp K) j , 

1 4^i 4fi 2 l pu J 


w: 


4^2 


L pu 


( 3 ) 

( 4 ) 


where = (1 + i) (up/2^ and j3 2 = (1 + i)(w/?/2// 2 )* , and K n are nth order modified 
Bessel functions and 9? denotes the real part of a complex number. The three constants A, B , C in 
(3) and (4) are given by the solution of 

AIo(jhR\) = Bio (/^Pi) + CKo (/?2#i) ? 

All (PiRi) = nllBhifaR^-CKrifaRi)], 

1 = BIq {P 2 R 2 ) + CKo (/?2^2) 5 


and 


A = i4[Ii(P2Ri)K o (02Ri) + Ki(02Ri)Io(P2Ri)] 


x 


| I\ (/3 2 Pi ) h (/?i Pi ) - Pi h (P2R1 ) A (^1 Pi ) 


Ko (P 2 R 2 ) 
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+ 


B = 


+ 


C = 


+ 


4 Jo (ftft) JJi (ftft) + 4 ft (ftft) Ji (ftft) Jo (ftft)} 

[4 Jo (ftft) ft (ftft) + 4 *0 (ftft) Ji (ft ft) 

{ [4 h (ftft ) Jo (ft ft ) - 4 Jo (ftft ) Ji (ft ft ) ft (ftft ) 

[4 Jo (ftft) ft (ftft) + 4 ft (ftft) Ji (ftft ) Jo (ftft)} 

[4 Ji (ft ft ) Jo (ft ft ) - 4 Jo (ftft) Ji (ftft) 

| [4 Ji (ftft ) Jo (ftft) - 4 Jo (ftft) Ji (ftft) ft (ftft) 
[4 Jo (ftft) ft (ftft) + 4ft (ftft) h (ftft) Jo (ftft) 


-l 


( 5 ) 


-l 


( 6 ) 


-1 


(7) 


We now non-dimensionalize fluid velocities by Wo, the basic steady flow evaluated at the axis of 
the core (see (3)) 


Wo = [r\ 0*j -V l) + JJImi] 


The length, time and pressure scalings are RiW 0 1 , and pW% respectively, and we define the 
following non-dimensional parameters: 


a = 


R P = 


#2 

WoR\p 

Mi 


m = 


A = 


M2 

Mi’ 

AJf^i J? e 


fl = 


ujR\ 

'Wo 9 


/J = (! + *) 


fLSl\ 2 


pWo 2 ’ 

The non-dimensional groups above represent the thickness of the undisturbed annular fluid, the 
viscosity ratio between film and core fluids, the non-dimensional frequency of the backgound oscilla- 
tory pressure gradient, a Reynolds number based on core variables, a non-dimensional amplitude of 
the backgound forcing and the Stokes layer thickness (proportional to |/3| — 2 . The non-dimensional 
form of the exact solution is 

/M_il e (inol 


Tft(M) = 1 - ™ r ‘ - 

a 1 4 - m — 1 
a 2 — r 2 


tft(M) = 


a 2 + m — 1 




+ » 

iA 


ttR 


J£L-[AI o (0r)-l]eW 

BIo(^t)+CKo 
\ m 2 / 


1 

m 2 


( 8 ) 

(9) 


In the expressions (8) and (9) the constants A,B,C are the non-dimensional versions of (5), (6) 
and (7) and are obtained by the replacements fi\ — > 1, fi 2 — * m -> P\R\ P and P2R1 — ► Py/ma and 

P 2 R 1 — 1 - y/rn(3. Note also that when A = 0 the steady problem is recovered as expected. 

Consider next a general unsteady axisymmetric CAF so that the velocity in regions j = 1,2 
is given by ([/, V, W)j (r, 2 , t) with the interface positioned at r = S (z,t). The Navier-Stokes 
equations in the core are 

» p. *-r p. n / rvO tt 1 Air o?. r r rr \ 

J , (10) 

( 11 ) 


dih 


rtr Ml 

, dPi 

1 

(d 2 lh 

1 du. 

a 2 ft 

dt 

rr 

dr 

+ W 'dz 

+ ~d 7 

Re 

^ dr ' 2 

r dr dz 2 

dWi 

rr dW l 

dWi 

, dPi 

1 

( d 2 Wj 

ldWi 

L d 2 w t 

dt 

+ ft - 

or 

+ W ' dz 

+ dr - 

Re 

l dr 2 

^ r dr 

dz 2 
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and in the annular region 


dU > + uS U ' 


dt 

dW 2 


dr 
dW 2 

+ u 2 —± + w 2 


dUn 

dP 2 

m 

(d 2 U 2 

ldU 2 

d 2 U 2 

1 dz ' dr 

R~ e 1 

\ dr 2 

r dr dz 2 

dW> 

dP 2 

m 

fd 2 W 2 

1 dW 2 

d 2 W 2 ' 

dz ^ dr 

~ Re 

\ dr 2 

^ r dr 

dz 2 


dt dr 

\ 

The equation of continuity in the core and annulus reads 

1 d(rUj) dW. 


dr 


+ 


dz 


= 0 , 


3 - 1 , 2 . 


( 12 ) 

(13) 

(14) 


The velocity (U 2 , V 2 , W 2 ) satisfies no-slip at the cylinder wall and (f/i, V \ , W\) must remain finite 
at the axis of the core. Using the notation 

[(•)]| = (*)i (O2 * 


R e j — 


Re j- 1 
R e m~ x j - 2 ’ 


the continuity of radial velocities, axial velocities, tangential and normal stresses at the interface 
are 


m 1 = 0, 

ml = 0, 


1 fdUi dWn 


Rej 

2 dU 


dz 


+ 


dr 


1 - 


05 V 

dz) 


+ 


2 (dUj dWjXds 


R, 


ej 


dr 


dz J dz 


= 0, 


R e j dr 




1 R e j dz 
J l d 2 S 1 


asy _2 

dz) + R, 


ej 


dUj dWj \ ds_ 

dz + dr J dz 


J 2 
T 1 


+ 


I ^-^1 1 + 


R 2 \ dz 2 


\dzj 


dS' 


j 2 

2\-i 


1+1 Tz) ) " °’ 


(15) 


where J = crpR\/ii\ is the non-dimensional surface tension parameter. The system is closed by 
imposition of a kinematic condition stating that the interface is a material surface, 


dS tt 7 dS 
Uj ~~at + w, lh 


at r = S (z,t ) . 


(16) 


3 The thin film limit 

The full nonlinear problem (10), (11), (12), (13), (15) and (16) presents a formidable analytical 
and computational task even under the reasonable assumption that the flow is axisymmetric. Con- 
siderable progress can be made, however, in situations when the annulus is thin compared with the 
core and when capillary instability is important (see Introduction for the physical significance of 
such limits.) We proceed, then, by formulating the stability problem for oscillatory core-annular 
flow where the annular layer is thin, referring to this layer as the film. In this regime a = 1 + 6 with 
f < 1 so that the film has thickness O (e) and we introduce a strained local coordinate y given by 

r— 1 + c(l — y). (17) 
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The flow described by (10)-(16) is now regarded as a perturbation to the basic state (U C) /, P c j), 
where the subscripts c and / denote the core and film regions respectively. The unknown interface 
position is written as r = S (z,t) = 1 + Srj (z,t), with 6 <C e (in the absence of a shear flow 6 
can be as large as 0(e) and the analysis presented below is still applicable) and to be determined 
by the main asymptotic balances (see below). The velocity perturbation in each layer is denoted 
by (U,V,W) c j and surface tension induces a jump in pressure P c j across the interface. The 
equation of normal stress indicates that this jump is of relative magnitude JSR~ 2 . The core and 
film dynamics are coupled by a balance of tangential stresses across the interface, and the relative 
asymptotic sizes of perturbed quantities in each layer then follow by balancing appropriate terms 
(see PMR also). It has been shown in PMR, and is also true for the present problem, that capillary 
forces affect the dynamics to leading order in the perturbation as long as the non-dimensional film 
depth, e, the surface tension group /, and the core speed characterized by the Reynolds number 
R e , satisfy the constraint eJ ~ R e . This is achieved by considering either of the following cases 


A: slow core flow with moderate surface tension, R e ~ e and J ~ 1; 

B: moderate core flow with large surface tension, R e ~ 1 and J ~ e -1 . 

To facilitate the balance of terms, a Taylor expansion about the unperturbed position r = 1 is 
used to describe the interface conditions. For example the undisturbed film velocity at the interface 
Wf(S(z, t)) may be written as the sum of the steady and unsteady parts 


W f (r,t) 
Wf (1 + Si)) 
W^ 8) (1 + Sr/,t) 


= W\ s) (r) + WY” (r,t) , 
26r] 


(us) 


2e £ 2 4e 

m m m 


m 


+ 


0(e 3 )+0( e «), 


= S 




mlo (/3) ml 0 ((3) 


1 + 


+ 


tn m-ai (3)) e20lH0) 

ph{p)) m/o(/3 ) 

0(fS)1 M: exp 


MW) 

min (0) 

+ 0(< 3 ) 


(«Ji 09) 


(ill£) j 


and the kinematic condition (16) becomes 

Oj = S^ + S (Wj + W,) g, at r = 1 + St). (18) 

It follows by leading order balances (see PMR for further details) that 

Uj ~ e 2 tf, 6 ~ £ 2 , r ~ Sty 


where r is the long time scale for the temporal instability. The oscillatory unsteady flow is of 
the form exp(il^), and in the analysis that follows we consider low frequency oscillations so that 
Sit ~ 1 with SI = € 2 (V Noting that 


P = 
h((3) = 
h((3) = 


(l + O*^)*, 
i + + o (Ry) , 


/ R e Slo \ 2 

[ iRSQo 0(^)1 

V 8 ) 

8 V / 


7 



we transform the coordinate system to one which oscillates and moves down the core with speed 
2em ~ 1 + O (c 2 ), by introducing the transformation 


- + --^ + o(A) 


2c 

rn 


€ 

m 


m i 


t - 


cA 


2 toO 


(1 4- c) + O sin (Oi) 


(19) 


With respect to this oscillating frame of reference the leading order kinematic condition reduces to 


_ 4 ~ d 7 ) ( 2 A ^ 


( 20 ) 


An evolution equation for rj(z,t) follows from (20) once Uf is expressed in terms of rj. This is 
achivable by matching with the core dynamics as described next. 


3.1 Slow flow, moderate surface tension 

Consider first the case when R e ~ e and J ~ 1. The perturbation velocities in the film and core 
become 


Uf 

— Uf — € 4 Ufo + • • • 5 

(21a) 

Wf 

= Wf + Wf = Wf + e 3 w f0 + ..., 

(21b) 

Pf 

= Pf + Pf = Pf + eP/o + • • • » 

(21c) 

Uc 

— U c = c 2 U c o + . . . , 

(22a) 

Wc 

= W c + w c = W c + € 2 W c0 + 

(22b) 

Pc 

= P c + P c = PcP tPcQ + • • • , 

(22c) 


and R e = Ae and A = 0 (1). Substitution of (21a-c) into (12)-(14) along with the local film scaling 
(17) yields the following leading order solutions in the film, 


Pfo 

W f0 

Ujo 


Pfo (2, t) , 

A |V d( f /o) 

m 2 dz 


+ yA(z 9 T) \ , 


A \y 3 9 2 (P J0 ) y 2 9 A 
m 6 dz 2 2 dz 


(23) 

(24) 

(25) 


The solution in the core is found in Fourier space and couples with the film through the balance in 
tangential stress to determine the function A (z, r). 


—m 


dWfp 

dy 


( 1 ) 


dU c0 

dz 


( 1 ) + 


dW t 


cO 


dr 


(i) 


(26) 


Substitution of (22a-c) into (10), (11) and (14) gives the following leading order core problem 


V 2 C/ c0 - ^ 

II 

^ 05 

o 

(27a) 

V 2 W c0 

,9P C o 
~ A dz ’ 

(27b) 

9(rUco) , 9W C o 

= 0. 

(27c) 


dr dz 
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which is to be solved subject to the continuity of velocity conditions 


U C 0 (r =1) = 0, 

" r -< r = 1 > = 2 "(^)- 

Defining a streamfunction if) in the usual way, U c o = — k W c o = f f^r, the system (27a-c) yields 
the familiar Stokes operator for if), 

( d 2 19 d 2 \ 2 , 

rdr + dz 2 ) ^ ~ °’ 


with solution readily obtainable in Fourier space in terms of modified Bessel functions 

^ r o° 

ip = ip(r,z)e~ tkz dz , 

J — OO 

= -rl 0 (k)F(k)rjI 1 (kr) + r 2 I 1 (k)F(k)rjI 0 (kr), 

F(k) = 


(28) 


2 (to — 1) 


TO [kit ( k ) - kI o (■ k ) + 2/o (A:) 7i (A;)] * 


The leading order radial velocity perturbation in the fi lm is then evaluated at y — 1 and on 
elimination of A(z,t) from (23)-(25), (26) and (28) can be written in the form 


2 i 


TO 


U f o(l) = N (k)r}+ —k'P f(h 


TO 


3to 


where the the kernel N ( k ) is given by 

N(k) = 


k 2 l\ (k) 


(29) 


(30) 


[kl 2 (k) - kl$ (fc) + 2kl 0 (*) h (*)] ' 

^From the kinematic condition (18) and continuity of normal stresses (15) respectively, we have 

U, 0 {y = i) = g-(£ + Acos(n„r)),g, (31) 


P/ ° = A2'( 7? + 


d 2 r\ 


and substituting (29) into (31) yields the desired evolution equation 


__l i [_ _ c os (SIqt) 7}— + r { ~o~o T ttt 1 + 


dr [ _(l _A_ 

dr \ to + 2 to 
i (to — 1 ) f+°° f+°° 


dr} J ( d 2 r} d 4 r/ 

dz 3toA i dz 2 dz 4 


7T7TO 


/ -l-oo ri-oo 

/ lV(fe)^(^ 1 ,r)exp(iA)(z — Z!))^ d/; = 0. 
-OO J — OO 


(32) 


It should be noted that (32) is a long-wave equation since it has been derived by assuming that 
the radial lengthscales in the film are much smaller than the axial ones. The validity of this 
assumption is checked a posteriori by studying the solutions of (32) and showing that no infinite 
slope singularities will arise. This has been proven for the KS equation analytically (see Introduction 
for references) and we provide evidence that this is also true for (32) from our extensive numerical 
work (see later). 
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3.2 Moderate flow, large surface tension 

The analysis presented in the previous Section applies also, with minor modifications when R e ~ 1 
and J ~ K For brevity we present the essential differences and refer the reader to PMR for a 
complete derivation. The dynamics in the film as well as the asymptotic solution there has the 
same form as above. Coupling between film and core is still achieved through the tangential stress 
balance and the normal stress balnace yields a pressure jump as before. The main difference is in the 
core dynamics: Since the Reynolds number is 0(1) now, the equations governing the perturbation 
there are the linearized steady Navier-Stokes equations (see PMR). The perturbation equations 
apear steady by virtue of the time-scale of the evolution. In particular, the core variables expand 
as 


U c — e 2 Uo -t- . . . , (33a) 

W e = W c + € 2 W 0 + ..., (33b) 

P c = 7 C + e 2 Po + (33c) 

The time-scale transforms according to the change of reference frame given by (19) and substitu- 

tion of this along with (33a-c) into the core Navier-Stokes equations (10)-(11) together with the 
intrduction of the streamfunction ip defined in Section 3.1, yields 


(i-r 2 )l z (m = 


D = 


R, 


-D 2 i>, 


®dr 2 - 
d r dr 



(34) 

(35) 


Solution of (34) is achieved by taking the fourier transform in z and the final evolution equation is 
obtained after satisfaction of all the boundary conditions (see PMR also). The result is 


dr] 



(& V flM 
y dz 2 dz 4 J 


i (m — 1) 
2k m 


/ +oo /*+ oo 

/ N (k)r](zi,T)exp(ik(z - z\))dz\ dk 

-oo J — OO 


0. 


(36) 


The kernel N(k) appearing above is expressible in terms of the confluent hypergeometric function 
and two auxiliary integrals given below, 


N(k) 

(k) 

N 2 (k) 


Ii(k)e~ x M(A, 2,2A) 
N^Ioik) - N 2 (k)h(kY 


f 

f 


(h(k)Ki(ks) - I 1 (ks)K 1 (k))s 2 e' Xs2 
(IoWK^ks) + h(ks)K 0 (k))s 2 e~ Xs2 


M(A, 2, 2\s 2 )ds, 
M(A,2,2\s 2 )ds, 


(37) 


where A = ^(kR e )^e~^ and A = 1 + - |, while M is the regular solution of the equation 


ZG” + (2 - 0& ~ AG = 0, 


(38) 


with primes denoting £- derivatives. Note that in order to solve (36) the kernel N(k) need only be 
computed once, given a value of v. This was done, in the absence of background fluctuations, by 
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PMR where it is shown numerically that the effect of the kernel (physically this is the viscosity 
stratification mechanism) is to organize the evolution into travelling waves for a wide range of 
parameters. 


4 Analysis of the evolution equations 


The time evolution of the interface position rj is described by equations (32) and (36). The following 
change of variables yields the canonical equation (39) given below 


r = 


ZmXti 


V — 


v 2 J 


- u , 


Jv 1 ■' 6A 

where L is the unsealed wavelength of the interface. 


z — v ix, v = 


L 2 ’ 


v^k\ 


du , .. du 

— -f (1 + cos(u;it))u— + 


d 2 u d 4 u 
dx 2 V dx 4 



(39) 


where 


Si 

CJl 

C{u) 


m 2 A 

4 ’ 
SmXilo 

Jv 
3 A 

TV Jv 


r r N (vh)u(x u t)e ikl -’ : -* 1 '>dx 1 dk. 

J — oo J —oo ' 


According to the change of variables above, we choose to construct solutions of equation (39) on 
spatially periodic domains; in physical space the limit v — > 0 corresponds to L — * oo. 

We begin to analyze the nonlinear dynamics of oscillatory CAFs by studying the oscillatory KS 
equation (referred to hereafter as the OSCKS equation). This limit corresponds to the nonlocal 
term in equation (39) being absent and is relevant if the viscosities of the two fluids are the same 
(m = 1) or if the Reynolds number is small for example of order e 2 or smaller (for a justification of 
the latter limit see [[8]]). The oscillatory KS equation (OSCKS) and the KS equation share the same 
linear stability spectrum when linearization is done about u = 0. More importantly f 0 * u(t,x)dx 
is a conserved quantity for both equations along with the fact that the energy equation for the 
evolution of / 0 27r u 2 dx is the same (this equation is easily obtained by multiplying (39) by u and 
integrating). As a result of this it can be concluded that when v > 1 solutions will tend to a 
uniform steady state and nontrivial dynamics first enter when v crosses below unity and linearly 
unstable modes are activated. 

The main purpose of this Section is to gain an understanding of the dynamics of the OSCKS 
equation. Due to the nonlinear nature of the problem we undertake extensive numerical work, but 
before presenting such results we consider the limit v — ► 1— which can be described analytically. 


4.1 Asymptotic solutions of the OSCKS equation for v near 1 

As mentioned above, when v > 1 there are no linearly unstable modes to pump energy into the 

system and if the initial condition has zero mean then the large time evolution yields a trivial steady 

state. ;From a normal mode linear stability analysis of the OSCKS equation (equation (39) with 
° 1 
m = 1), it is clear that that the number of linearly unstable modes is mod [v~i] and by setting 

V = 1 - 6i, 0 < €i « 1, 
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we are in a position to study the evolution in the neighborhood of the first unstable mode. In 
what follows we begin by consideration of solutions with zero mean and odd parity. We note that 
for v near 1 (but not necessarily asymptotically close), the KS has a global steady state attractor 
which can be cast into an odd-parity zero-mean profile by Galilean transformation and translation 
invariance. Using Fourier analysis, then, we can represent zero mean odd-parity solutions by the 
infinite sum 

oo 

u(x,t) = ^2 ak(t) sin (kx). (40) 

k~l 

The partial differential equation is dissipative and we expect a low modal behavior as for the KS 
equation. In practice this means that a Galerkin approximation can be employed to truncate the 
sum in (40) to N terms, with N depending on v and also on tfi. £From previous work on the KS 
(see PS and SP) we find that a few modes more than the number of linearly unstable ones provide 
a sufficient description of the attractors. In the analysis that follows the asymptotic limit €i — ► 0 
is described by the two leading modes a\ and a 2 . The equations satisfied by a\ and are 

(41a) 
(41b) 

This is a two-dimensional dynamical system and our main concern is with the behavior at large 
times which provides the features of the attractor for this range of values of v. In the context of 
the Galerkin approximation the system (41a-b) is exact; its validity hinges on the fact that there 
are only two modes of importance, something which can be true if v is near 1 and €\ is small. 
In practice ci need not be asymptotically small for the two-mode truncation to be a reasonable 
approximation; this is shown later (see Figure 2) by a computation with — 0.1 and with 2 modes 
or 20 modes respectively with indistinguishable results. The smallness of can be used, however, 
in an analytical description of these solutions and consequently provide an accurate description of 
the dynamics. 

In order to motivate the multiple scales ansatz to follow, we consider first the solution of (41a-b) 
when = 0. If the unsteady terms are dropped the steady states are obtained by solution of two 
coupled algebraic equations which yield a\ = y/48ei and a 2 = —2e\. There does not appear to be 
a closed form solution for a\ and 02 at general times but a multiple scales analysis can be used to 
provide the essential features of the dynamics. It can be seen that since a\ ~ yJT{ and 0-2 ~ «i, the 
principal time scale is long and of order Cj 1 . With <5i non-zero (and of order 1 or smaller) the 
appropriate expansion is 

aj(J) = c} /2 [Aio(ti,t 2 , . . •) + €ij4n(/i,t2> •••) + •••], (42a) 

a 2 (0 = «l [i42o(*l»*2» • • •) + €iA2l($l,t2r •••) + •••]» (42b) 


da i , f , .\\ a i a 2 n 

— - - (1 + tfi cos(wit))— eiai = 0, 

at £ 


Ip- + (1 + #1 cos(w,«))^ + (12 - 16 c,)a 2 = 0. 

at 2 


where the multiple scales are given by 

t\ — t, ^2 — ^3 — • • • • 


Substitution of (42a-b) into (41a-b) yields the following equations to leading order 


d^4io 

dti 


= 0, 


(43) 
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( 44 ) 


- -(1 4- cos(u>ifi))AioA 2 o - Aio = 0, 
q 4- -(1 + cos(u;iii))A 2 0 + I2A20 = 0, (45) 

-777— 4 — tt~ 4- (1 + ^1 cos(u>i£i))Aioj 4 h + I2A21 — I6A20 = 0* (46) 

Ot\ ot 2 

It can be seen from (43) that A\q is independent of t\ and it follows from (45) that 

A20 = K\e 12<1 — A\ 0 ^67 cos(wiii) + - 7W1 sin(u;ifi) + , (47) 


where 


h 

144 4- ' 


The value of K\ in (47) is fixed by the initial conditions and does not affect the large time behavior 
of the solutions. With .A20 available, then, we turn to (44). Since Am is independent of fi, solutions 
of (44) will become unbounded as t\ increases unless a secularity condition is satisfied. This arises 
by setting to zero terms in the forcing function which do not involve £1; this is shown explicitly 
below by a re-grouping of terms in (44). 


- \a w Ke~ 12t '( 1 + S 1 cos(wjti)) - (67 + ^)A? 0 cos(w 1 < 1 ) 

1 1 
-^17^10 sin^fi) - 37<5iA 2 0 cos(2wifi) - -u\i8\A\ 0 sin(2wi$i) 

+ “ 2a '° ~ Al ° = °’ 


where 


2 _ 1 1 3 7 ^i 

" “ ' ,48 + ~ 


q;^Ajo — A 10 — 0, 

for Am, 

A 2 - ^ 

10 a 2 L + e~ 2f 2 ’ 

where K 2 is fixed by the initial conditions. With Am and A 2 o known, higher order terms can be 
calculated by direct integration and application of secularity conditions where necessary. 

Of particular interest is the comparison between the asymptotic solutions just described and 
the numerical solutions of (41a-b) for small C \ . For large times the leading order solutions are 


A 10 ~ 

(s 

+ 2(144 + w?)J ’ 

(48a) 

An ~ 

1 

2a 3 

(j + 24a, 1 ) Sm( “ ltl)+ 2 C ° S( “ ltl) 

(48b) 

- 

37^>i 

2c^i 

7^1 

sin ( 2 a;i 1 1 ) -I — — cos( 2u>i ti ) , 
8 

(49b) 


a 20 ~ -(js + apJtw?) ) . (49c) 


The secularity condition is, then, 

dAio 
dt 2 

which gives the following explicit solution 


+ 
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It is seen from (48a-c) that the leading order solution Aio of a i is a constant while the next 
order correction, namely An, provides an oscillatory component of period 2rc /uj\. For the Fourier 
mode a 2 the oscillations are captured at leading order by the asymptotic expression A 2 o the period 
being 2ir /u\. These analytical results suggest, therefore, that the effect of the oscillatory pressure 
gradient in the undisturbed flow is to induce a nonlinear interfacial evolution (within the framework 
of the present work) which is (i) spatially periodic with wavelength 27 rz /~2 with v slightly less than 
unity, and, (ii) which is time-periodic with frequency locked to that of the driving pressure gradient. 
Such behavior is also seen in related two-dimensional systems such as the forced Duffing oscillator 
for instance (see Chow h Hale (1982) for example). The forcing in the present problem is different 
in that it only alters the distribution of energy rather than adding or removing energy from the 
system. The evolution of the OSCKS becomes increasingly complicated as v decreases and the 
dimension of the attractors increases; in the following Section we carry out extensive numerical 
experiments in order to gain an understanding of the attractors for this system. 

Before presenting direct numerical simulations of the OSCKS, we provide a comparison be- 
tween the multiple scales asymptotic theory and the numerical solution of (41a-b). The important 
parameter is and numerical results with different values of uq and 8\ are qualitatively similar; 
in what follows, therefore, we fix the values uq = 1 and £1 = 0.5 and carry out an evaluation of 
the asymptotic theory for a range of t\. Equations (41a-b) were integrated numerically by spec- 
ifying initial conditions eq = 1, a 2 = 0 at time t — 0 and solving to large enough times so that 
any transient behavior disappears. The Fourier component a\ is compared with the corresponding 

1 3 

two-term asymptotic solution ef A 10 + ef A u , while the second Fourier component a 2 is compared 
with its corresponding one-term asymptotic solution €\ A 2 o- Numerical solutions are given in Figure 
1 for values of — 0.1 and = 0.01. The agreement is seen to be very good with additional 
improvement as €\ decreases. By analyzing the phase plane of the energy, we have verified that 
the solutions are indeed periodic and frequency-locked to the background flow oscillations, with a 
period of 2ir for both cases. The fact that the two-mode truncation is a very good approximation to 
the solution is illustrated in Figure 2. A numerical solution has been obtained using 20 modes (see 
the next section for details) and a comparison is made between the Fourier coefficients a 1 and a 2 as 
well as the corresponding asymptotic expressions. The dotted lines correspond to the asymptotic 
results and it is seen that the 2 mode and 20 mode truncations produce almost identical results. 

5 Numerical experiments on the OSCKS 

The solution of the full nonlinear evolution equation (39) for general values of v , 8\ and aq must 
in general be obtained numerically. The method we apply here is an adaptation of the scheme 
developed by PS and SP for the solution of the KS equation, which is obtainable from (39) when 
m = 1 and £1 = 0. 

5.1 Numerical methods 

The following Galerkin approximation is used to approximate solutions with zero mean 

N 

u{t,x) = ^2 ( a k(t) sin (kx) + bk(t) cos (kx )) . (50) 

k = 1 

When m / 1 the solution contains both sines and cosines whereas when m — 1 the symmetry of 
(39) allows odd-parity solutions for which 6* = 0 in (50). Such odd-parity solutions would evolve 
from odd-parity initial data for instance; general initial conditions, however, require both sines and 
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cosines and so there are twice as many unknowns then. Briefly, then, when (50) is substituted into 
(39) and coefficients of the different harmonics are set to zero, a nonlinear system of 2 N (or N) 
coupled first order differential equations are obtained for the evolution of the fourier coefficients 
(see PS). These equations when written as a first order system, contain a nonlinear coupling term 
and a linear term resulting from the fourier projection of the linear spatial operators in (39). A 
split-step scheme is implemented with the linear part (which contains all the stiffness at large N and 
the dispersive effects due to the pseudo- differential operator) integrated exactly and the nonlinear 
part with a fourth order variable step Runge-Kutta method. Accuracy tests are applied at each 
time-step. For particular experiments the number N is chosen large enough so that the minimum 
amplitude of any of the fourier components is below a tolerance level (this tolerance is usually set 
to < 10 -8 ). The linear dispersive term also imposes a restriction on the time-step as we explain 
next with the model ut -f- u x xx — 0* The solution of the kth. fourier coefficient after a time-step 
At is proportional to exp(*fc 3 At); the exponential is evaluated accurately, then, if k^^At <C 1- In 
practice we allowed the upper bound to be about 0.1 — 0.3 without loss of accuracy. All numerical 
experiments described below and which contain high frequency oscillations superimposed onto the 
lower frequency background flow have been confirmed by convergence studies in both N as well as 
the time- step. 

In the presence of periodic forcing (<!q ^ 0) solutions of equation (39) display a variety of 
interesting behavior as v decreases below 1. For various values of uq and <$i we have characterized 
the effect of the oscillations on a number of attractors. In Table 1 we provide some representative 
behavior for odd-parity solutions generated by the initial condition w(0,£) = — sin(&); these results 
are discussed in Section 5.2 below. General initial conditions were also used to compute asymmetric 
profiles; in many instances the initial conditions are taken to be random and so the attractors we 
compute and describe below in Section 5.3 and Table 2 are those with the largest basin of attraction. 
In Section 5.4 we give some results for m ^ 1 and in particular try to quantify the effect of large 
dispersion on the dynamics. We use various diagnostic tools to analyze the data of our numerical 
experiments and use them to infer the type of attractors present. These techniques are summarized 
in the Appendix. 

5.2 Odd parity solutions of the OSCKS 

When = 0, odd-parity solutions of the KS equation for 1 > v > 0.06 (approximately) are 
attracted to a non-uniform steady state. Within this window, are sub- windows where fully modal, 
bimodal and trimodal attractors are observed, with respective spatial periods of 27 r, 7 r, and 37r/2. 
For ^ / 0 and for 1 > v > .06 our numerical experiments show that solutions of the OSCKS 
equation lock into the frequency of the forcing pressure gradient S x cos (uqi). 

For example, the attractor of the odd-parity KS at v — 0.1 yields a bi-modal steady state which 
becomes time-periodic when ^ 0. The time dependence is seen in Figure 3a where the evolution 
of the energy of the solution, E(t) = £ / 0 27r u 2 (t, x)dx say, is plotted for two different frequencies, 
aq = 1,2. It is found that E{t) is periodic of period ^ . The spatial characteristics of the flow 
are preserved by the forcing as shown Figure 3b. This picture was produced as follows: If the flow 
is periodic in t, then u(t,x ) = u(t + j^,®) and so plotting u(t,x) at time intervals separated by 
22L should yield the same curve as seen in Figure 3b. Further more we chose to plot u at times in 
its cycle where E(t) is maximum and u x (t, 0) > 0 and compared the profile with the steady-state 
solution of the KS for 6\ — 0, finding identical results. It can be concluded, therefore, that the 
forcing in this case produces a standing wave locked onto the driving pressure gradient. 

The KS equation also supports windows in the parameter space ( v ) where chaos is preceded by 
a complete sequence of period doubling bifurcations (see SP and PS). One particular example is 
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the case v = 0.0303 for which the period is 0.89 time units. The behavior of this attractor in the 
presence of harmonic background forcing was analyzed by following the evolution of u for various 
values of and u > Even for small <?q the periodic attractor is found to become quasi-periodic 
after a long time. Figures 4(a-c) show E(£), the phase plane of E(t ) and the return map of the 
minima, for increasing £1 = 0.01,0.1,0.5 respectively, and fixed forcing frequency oq = 1.0. As <!q 
increases from 0.01 to 0.1, it is seen from Figures 4(a) and 4(b) that the energy has two dominant 
time scales: an oscillation with period roughly that of the KS with = 0 and a longer modulation 
with period roughly 2ir provided by the forcing pressure gradient. The interaction between the two 
frequencies produces a quasi-periodic flow as can be seen from the return maps in Figures 4(a), (b) 
which are two closed unfolded loops in both cases. One of the loops in Figure 4(a) is very flat and 
appears as a line but an enlargement shows an open loop. The phase-plane is more revealing for the 
<$i = 0.01 case; the lower modulational frequency is responsible for the larger turns while the higher 
frequency produces the smaller turns in the phase-plane. The spacings between turns eventually 
fill out more and more indicating the quasi-periodicity confirmed numerically by the return maps. 
When <!q is increased further to a value of 0.5, it is found that after approximately 160 time units 
of irregular modulated high frequency oscillations the solution locks into the forcing frequency and 
its period is 27 r. This is confirmed further by the phase-plane and the return map which contains 
only a single point - see Appendix. 

For v <C 1 we observe regions of chaotic oscillations for the solution of both the KS equation and 
the OSCKS equation. The aperiodic oscillations alternate with periodic/ quasi-periodic regions, but 
the lengths of the subwindows become vanishingly small as v decreases. A phase-plane analysis and 
the return map of the minima do not indicate recognizable behaviour, such as strange attractors. 

5.3 Solutions of the OSCKS equation with general initial data 

We have also computed a number of solutions of (39) for general initial data. These computations 
are more expensive that those of Section 5.2 since there are twice as many equations now. This class 
of solutions gives richer dynamics for both the KS and OSCKS equations. For the KS equation, 
for example, we observe various steady state attractors for a wide range of v < 1, but in addition, 
additional dynamics are found such as travelling waves and homoclinic bursting (periodic and 
aperiodic). Consider for example the case v = 0.3 for which the KS equation yields a travelling wave 
solution. Figure 5(a) shows that the energy quickly evolves to a constant value (4.27 approximately) 
and a steady-state travelling wave emerges after a long time whose phase-speed is small but non- 
zero. The travelling wave can be seen from the last figure which was produced by recording the 
solution after 11 equal time intervals and shifting each profile verically by 5 units. Figures 5(b), (c) 
show the effect of the background oscillations on the steady-state travelling wave depicted in Figure 
5(a). The values of 6\ are 0.01 and 0.1 respectively and the OSCKS equation is found to yield a 
wave which oscillates quasi-periodically and also travels. The return maps indicate that the solution 
is quasi-periodic. 

In the window 0.23249 > v > 0.17735 solutions of the KS equation undergo periodic homoclinic 
bursting. A representative case has v = 0.22 and is depicted in Figure 6(a) which shows that 
the energy remains constant for approximately 21 time units between each localized burst. The 
phenomenon is time periodic and repeats itself as indicated by the phase plane which contains 
about 20 turns (for the developed solution) superimposed on each other. Even for a small 6 = 0.05, 
the solution of OSCKS equation appears to be chaotic even after a significant time has elapsed. 
The phase plane analysis and the return map do not give convincing evidence of quasi-periodicity. 
It is interesting to note from the return map that the solution is almost quasi-periodic but drifts 
away from the what would have been closed loops. 


16 



For smaller solutions of both the KS and OSCKS equations show mostly chaotic oscillations. 
We present a typical case for the KS with v = 0.1212; the time history of E(t) is chaotic but the 
phase plane is not sufficient to determine the type of oscillations present. The return map in Figure 
7(a), however, indicates the presence of a strange attractor made up of live branches joined in a 
common region where folding takes place. Figure 7(aa) presents graphical evidence that the strange 
attractor is self-similar, by carrying out successive enlargements of the original return map; the 
dimensions of the first picture are 0.12 x 0.12 and by the fifth enlargement the size is 0.0001 x 0.0001 
so that an enlargement of 1200 has taken place with the self- similarity clearly apparent. In the 
results presented in Figure 7(b) we quantify the effect of a background oscillation with Si = 0.1 
on the chaotic motion of Figure 7(a). It can be seen that E(t ). is highly irregular and its phase 
plane yields a picture confirming this. The return map fails to produce a recognizable pattern; in 
fact from the return map we can exclude the possibilities of either quasi-periodic motion or chaotic 
motion with a strange attractor similar to the 0 solution. It is possible that a strange attractor 
in a higher dimension than the plane is present (such a possiblity could be studied by forming the 
return map of triplets etc.) but this is not explored here. 

5.4 Effects of viscosity stratification, (m^l) 

When the viscosities of the film and core fluids differ, a non local term L(u) is introduced into 
the evolution equation (see (39). In the absence of oscillatory modulations in the driving pressure 
gradient, it has been shown in PMR that even for relatively small viscosity stratification (i.e. for 
m near 1), otherwise chaotic solutions of the equation describing the interfacial position (given by 
numerical solution of the KS equation), are organized into regular travelling wave pulses. This 
is true for both the canonical regimes presented in Sections 3.1 and 3.2 respectively. Physically 
this reflects the effect of dispersion provided by the viscosity stratification on the nonlinear waves; 
it has been established numerically for a wide range of v that for = 0 there exists a critical 
value of m (which depends on v) above which travelling waves are attained. The capillary terms 
(second and fourth derivatives) act to set the scale of the pulses which have been found to increase 
inamplitude as the dispersion increases. A detailed study of dispersive effects in core- annular flows 
can be found in Coward, Papageorgiou k Smyrlis (1994). To illustrate these results we present a 
numerical solution with the piirely disspersive Bessel functions kernel of Section 3.1 at a value of 
v = 0.1212 and = 0 which for m = 1 yields chaotic oscillations (see Figure 7a). For m = 3.0 the 
large time evolution is attracted to a travelling wave state as shown in Figure 8 which depicts the 
time evolution over 4.5 spatial periods, by plotting the profiles at successive equispaced times by 
shifting them vertically by a small amount to facilitate the graphical representation. 

A forced oscillation is now imposed on the flow of Figure 8 with a 8 = 0.1. The results are given 
in Figure 9: From the energy plot we see that a small oscillation persists about a constant value. 
This oscillation is brought out more clearly in the phase plane plot which indicates high frequency 
oscillations modulated by a low frequency motion provided by the forcing. A convergence study 
was done to verify that the high frequency oscillations are present and not due to numerical errors. 
In fact the modulated high frequency response is shown to be quasi-periodic by the return map 
of the energy minima shown in the last of Figure 9. The nonlinear travelling wave in the absence 
of forcing, then, is found to still preserve its travelling character as expected but also oscillates in 
time quasi-periodically. A picture of the wave motion is provided in Figure 10 which was produced 
in the same way as Figure 8 with a vertical shift of 3 units. The travelling component of the wave 
is easily discerned; the quasi-periodic modulations cause the small displacements in the wave and 
consequently the dark lines which are an indication of the locus of the wave minima in the x — u 
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plane are not straight lines but modulated around a straight line. An observation which is clearly 
illustrated by the spatial solution of Figure 10 in which u is given at successive time intervals, with 


6 Conclusions 

A weakly nonlinear asymptotic theory has been used to derive some novel partial differential equa- 
tions that govern the stability of two-phase core-annular flow when an oscillatory pressure gradient 
acts. The equations extend the familiar Kuramoto-Sivashinsky models of CAFs to analogous evo- 
lution systems with time periodic coefficients. The main interest of this work is in the study of 
interfacial instability of oscillatory two-phase flows and the parametric evaluation of its effects on 
conventional CAF caused by constant pressure gradients. Such studies are desirable in the deter- 
mination of physical mechanisms in complicated fluid flows which may be used advantageously to 
control inherent instabilities. 

The main findings of our numerical experiments can be summarized in the following categories: 
(i) Stationary steady-states become time-periodic with frequencies locked to that of the forcing; (ii) 
Time periodic states become quasi-periodic, chaotic or lock into the forcing frequency if the forcing 
amplitude is large enough (for odd-parity solutions at least), (iii) Travelling wave states bacome 
travelling quasi-periodic waves, (iv) Periodic homoclinic bursts become unrecognized chaotic states, 
(v) Chaotic states which possess self-similar strange attractors become unrecognized chaotic states. 

The forcing appears to promote irregular oscillations in regimes with large dispersion which 
would otherwise produce travelling wave pulses in unforced flows. The solutions in this case are 
perturbed by a quasi-periodic in time modulation whose amplitude increases with the forcing 
amplitude. 

Our numerical results have not revealed any situations where the solutions lock into subharmonic 
or other lower frequencies with respect to the background forcing oscillation. It appears that if 
there is frequency locking this is at the driving frequency, for a wide range of parameters. 
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A Analysis of the numerical data 

The diagnostics we use are aimed at the determination of the type of attractors present. To do this 
we first begin with the L 2 -norm (or energy) of the solution define d by 

1 /‘ 27r 

E(t) = — u 2 (t,x)dx , 

I'K Jo 

and which can be computed along with the solution. For example, in the odd-parit y case we have 

z 1 

with analogous expressions in the general case. Of concern is the behavior of E(t) for large t. If 
E(t) tends to a non-zer o constant as t — * oo, for example, then u(t,x ) tends to a non-uniform 
stationary or travell ing steady state. Time periodicity in the fourier coefficients (and consequently 
in u(t,x)) show s up as a time periodicity in Eft). A useful way to confirm such periodic flows is 
to constru ct numerically the phase-plane of (E,E); the index of the closed curve determines the 
numbe r of maxima and minima of E(t) over one period and can be used with good accuracy to 
follow pe riod-doubling bifurcations, for instance, as in (SP) where as many as thirteen bifurcations 
we re analyzed this way. 

When the dynamics gets more complicated and chaotic, the phase plane may gain mo re and 
more turns as the dynamics enter a strange attractor, fill out whole regions of the (E,tE) plane 
as t increases as for quasi-periodic oscillations for instance, or produce dynami cs which are both 
random and unrecognizable. In all these cases it is difficult to infer with much certainty what is 
happening and a more accurate test is needed. One diagnostic we found to be very useful is the 
return map described below. 

A return map can be constructed as follows. Given the time history of E(t) (th e non-constant 
state is the interesting one), and the phase plane described above, Poincare secti ons can be 
constructed in the usual way. In particular the minima and maxima of E(t) are Poincare sections 
having E = 0 and negative and positive respectively. For def initeness consider the minima of 
Eft) from a given point in time (large enough for transients to be unimportant) onwards and define 
them by the infinite sequence {J5 n }^°. The return map of these values, which gives a picture of the 
attractor in that i t is a construction of the map (it may not necessarily be a function) necessary to 
reproduce the flo w, is obtained by collectively plotting the pairs points {( Ek , E k+ i)} k=1 fty in the 
plane. In practice we do not compute an infinite number of minima but must nonetheless compute 
enough to provide a picture of the dynamics. In addition is is essential that the values o f the 
minima be found with enough accuracy from the computation - see accuracy requirements belo w. 

The simple example of the logistic map serves as a useful illustration of the ab ove construction. 
The logistic map produces a sequence of iterates {z n }^° from the quadratic formula x n +i = 4fj,x n {l — 
x n ) = f(x n ) with 0 < n < 1 a parameter. A gr aphical iteration starting from some initial value 
produces the desired sequence (see Devaney for instance) . This iteration and proceeds as follows: 
Start from a given point x k say; the iterate of x k is x k+1 = f{x k ) whose value is the length of the 
vertical line joining x k to the curve; the iterate of x k+ i is found by moving horizontally from the 
point (x k ,x k+1 ) to the line of unit slope through the origin, and then vertically up to the curve 
and the process can be repeated ad infinitum. Collectively, then, the points {(£&> x k +\)} k=l are 
drawn in the plane and all lie on the parabola of the logistic map. The construction of the return 
map from a sequence such as duplicates exactly this graphical iteration and constructs 

the unknown curve f(x), if it exists, or a more complicated object such as a strange attractor for 


21 



instance. Snch return map diagnostics have been used extensively by experimentalists - see Berge 
et al. for numerous examples. 

The interpretation of the return map is usually quite clear provided the minima have been 
determined accurately. There are five distinct pictures which we have encountered in our numerical 
work on the KS and the OSCKS based on the r eturn map of the minima of E(t) and which are 
interpreted as follows as far as the dynamics of E{t) are concerned: 

1. One or more fixed points in the plane - time periodic flow with a number o f minima equal 
to the number of fixed points. 

2. An open curve in the plane which may or may not be a single valued functio n - chaotic 
flow as would be obtained from the logistic map for instance (the former case), or ju st 
after the Feigenbaum cascade to chaos for the KS equation (see PS, SP). 

3. An open curve which is highly folded and which exhibits self-similar behav ior at smaller 
scales - chaotic flow with a strange attractor of dimension between 1 and 2 reme niscent 
of the Henon map - such behavior has been found by Papageorgiou k Smyrlis (1994) f or 
the KS equation far from the Feigenbaum cascade. 

4. One or more closed curves with no folding - quasi-periodic flow as reporte d in several 
examples in the present article. 

5. Unrecognizable pattern which fills out regions of the plane - we term this type of flow as 
unrecognizable chaos and several examples are given here and in Papageorgiou k Smyrlis 
(1994) for the KS equation. 

Finally we make a comment about the need of high accuracy requirements. As an il lustration of 
such needs, consider a case which produces quasi-periodic oscillations in E(t). E(t) is numerically 
computed at intervals separated by the time-step At. In many problems a large At can be used to 
compute the flow accurately. The simplest approach of finding the energy minima is to dynamically 
track the minimum value and store it when found; the error associ ated with this is of the order of 
magnitude of the local value of E x At] instea d of a closed curve of zero thickness such errors can 
indicate spurious behavior which can be of the same size as real phenomena (this is especially true 
as becomes increasingly small er and the amplitude of quasi-periodic modulations decreses). 
To overcome this difficulty the minima are computed dynamically in the code by a high order 
interpolation involving as many as 15 points, depending on the time-step, on either side of a local 
maximum or minim um . All results reported here have been produced by an accurate evaluation of 
the minima and maxima of E{t) in order to produce highly accurat e and convincing return maps. 
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Tables 


V 

N 

h 

U>1 

Description 

Figure 

0.9 

10 

0.0 

0.0 

Uni-modal steady attractor 


0.9 

10 

0.5 

1.0 

Uni-modal periodic attractor, T — 2tt 

2 

0.1 

16 

0.1 

1.0 

Bi-modal periodic attractor, T = 2tt 

3 

0.1 

16 

0.1 

2.0 

Bi-modal periodic attractor, T = ir 

3 

0.1 

16 

0.1 

0.1 

Bi-modal periodic attractor, T = 7r / 5 


0.0303 

20 

0.0 

0.0 

Periodic attractor 


0.0303 

20 

0.01 

1.0 

Quasi-periodic oscillations 

4(a) 

0.0303 

20 

0.1 

1.0 

Quasi-periodic oscillations 

4(b) 

0.0303 

20 

0.5 

1.0 

Periodic attractor, T — 2ir 

4(c) 

0.0215 

24 

0.0 

0.0 

Chaotic oscillations 


0.0215 

24 

0.1 

1.0 

Chaotic oscillations 


0.0215 

24 

0.5 

1.0 

Chaotic oscillations 



Table 1 : Overview of solutions of the OSCKS with sinusoidal profile. 


V 

N 

Si 

U>i 

Description 

Figure 

0.3 

12 

0.0 

0.0 

Traveling wave 

5(a) 

0.3 

12 

0.01 

1.0 

Quasi-periodic traveling wave 

5(b) 

0.3 

12 

0.1 

1.0 

Quasi-periodic traveling wave 

5(c) 

0.3 

12 

0.5 

1.0 

Quasi-periodic traveling wave 


0.22 

16 

0.0 

0.0 

Periodic homoclinic bursting 

6(a) 

0.22 

16 

0.1 

1.0 

Chaotic oscillations 

6(b) 

0.136 

16 

0.0 

0.0 

Bi-modal steady attractor 


0.136 

16 

0.0 

1.0 

Modulated bi-modal attractor 


0.1212 

20 

0.0 

0.0 

Chaotic oscillations/Strange attractor 

7(a) & 7(aa) 

0.1212 

16 

0.1 

1.0 

Chaotic oscillations 

7(b) 

0.1212 

16 

0.5 

1.0 

y Chaotic oscillations 


0.085 

16 

0.0 

0.0 

Chaotic oscillations 


0.085 

16 

0.1 

1.0 

Chaotic oscillations 



Table 2 : Overview of solutions of the OSCKS with mixed profile. 
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gure ' 


00 15 


epsilon=0.1 


epsilon=0.01 


Comparison of asymptotic, and numerical solutions of the OSCKS equation with v near 
1 ? = 0.5, uj\ — 1.0 . The fourier coefficients a\{t) and a 2 (t) for e = 1 — v — 0.1,0.01; 

dark line - asymptotic result (48a-c). 













3 OSCKS equation: Bi-modal periodic attractor, with v = 0.1, S 1 - 0.1, - 1, 

Energy against time, and the spatial structure of the profile. 
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4(c) OSCKS equation: quasi-periodic oscillations initially and later frequency locking to 
time-periodic state, with v — 0.0303, £1 = 0.5, = 1.0. The KS in 4(a)- (c) has a 

time-periodic attractor of period ~ 0.89. 









Figure 5(a) 






Energy 


Figure 5(b) 

4.4i 1 1 1 i T 



Return Map of Minima 



5(b) OSCKS equation: quasi-periodic travelling wave, v — G.3, <?>i — 0.01, — 1.0. 


































7(a) KS equation: Chaos just beyond a period-doubling cascade; strange attractor, 
0 . 1212 . 





7(aa) Enlargement of returm map sh< 
v — 0.1212. Overall enlargement 
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Energy 



Figure7(b) 
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Figure9 



9 Modified OSCKS equation with viscosity stratification, m = 3.0, v = 0.1212, S\ — 0.1, 
cji = 1.0. Quasi-periodic attra ctor. 
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Figure 10 


200 h 


150h 



lOOh 


10 Modified OSCKS equation with viscosity stratification: spatial evolution of quasi- 
periodic travelling wave, m — 3.0, v = 0.1212, Si — 0.1, uji = 1.0. 
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